mod2a <- lm(sadness_any ~ t1 + block, data = dat2, weights = ipw)
mod2a$se <- coeftest(mod2a, vcov = vcovHC, type = "HC1")[,2]
mod2b <- lm(sadness_any ~ t1 + block + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod2b$se <- coeftest(mod2b, vcov = vcovHC, type = "HC1")[,2]
mod3a <- lm(fear_any ~ t1 + block, data = dat2, weights = ipw)
mod3a$se <- coeftest(mod3a, vcov = vcovHC, type = "HC1")[,2]
mod3b <- lm(fear_any ~ t1 + block + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod3b$se <- coeftest(mod3b, vcov = vcovHC, type = "HC1")[,2]
stargazer(mod0a, mod0b, mod1a, mod1b, mod2a, mod2b, mod3a, mod3b,
se = list(mod0a$se, mod0b$se, mod1a$se, mod1b$se, mod2a$se, mod2b$se, mod3a$se, mod3b$se),
digits = 2, no.space = T, keep.stat = c('n', 'rsq'),
covariate.labels = c('Anger Appeal', 'Members (Log)',
'Any Violence', 'Poverty',
'Constant'),
omit = c('block'))
#####
## Table 2 - outcomes
#####
## ITTs
dat2$pre <- dat2$index.pre.m
mod0a <- lm(index.m ~ t1 + block, data = dat2, weights = ipw)
mod0a$se <- coeftest(mod0a, vcov = vcovHC, type = "HC1")[,2]
mod0b <- lm(index.m ~ t1 + block + pre + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod0b$se <- coeftest(mod0b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$msgs.pre
mod1a <- lm(msgs ~ t1 + block, data = dat2, weights = ipw)
mod1a$se <- coeftest(mod1a, vcov = vcovHC, type = "HC1")[,2]
mod1b <- lm(msgs ~ t1 + block + pre + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod1b$se <- coeftest(mod1b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$resp.pre
mod2a <- lm(resp ~ t1 + block, data = dat2, weights = ipw)
mod2a$se <- coeftest(mod2a, vcov = vcovHC, type = "HC1")[,2]
mod2b <- lm(resp ~ t1 + block + pre + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod2b$se <- coeftest(mod2b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$viva.pre
mod3a <- lm(viva ~ t1 + block, data = dat2, weights = ipw)
mod3a$se <- coeftest(mod3a, vcov = vcovHC, type = "HC1")[,2]
mod3b <- lm(viva ~ t1 + block + pre + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod3b$se <- coeftest(mod3b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$vic.pre
mod4a <- lm(vic ~ t1 + block, data = dat2, weights = ipw)
mod4a$se <- coeftest(mod4a, vcov = vcovHC, type = "HC1")[,2]
mod4b <- lm(vic ~ t1 + block + pre + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod4b$se <- coeftest(mod4b, vcov = vcovHC, type = "HC1")[,2]
stargazer(mod0a, mod0b, mod1a, mod1b, mod2a, mod2b, mod3a, mod3b, mod4a, mod4b,
se = list(mod0a$se, mod0b$se, mod1a$se, mod1b$se, mod2a$se, mod2b$se, mod3a$se, mod3b$se, mod4a$se, mod4b$se),
digits = 2, no.space = T, keep.stat = c('n', 'rsq'),
covariate.labels = c('Anger Appeal', 'Pre-Treatment Outcome',
'Members (Log)',
'Any Violence', 'Poverty',
'Constant'),
omit = c('block'))
#####
## Text Section 5.2 - percent change in treated
#####
ates <- c(mod1b$coef['t1'],mod2b$coef['t1'],mod3b$coef['t1'],mod4b$coef['t1'])
means1 <- c(mean(dat2$msgs[dat2$t1==1],na.rm=T), mean(dat2$resp[dat2$t1==1],na.rm=T),
mean(dat2$viva[dat2$t1==1],na.rm=T), mean(dat2$vic[dat2$t1==1],na.rm=T))
means0 <- c(mean(dat2$msgs[dat2$t1==0],na.rm=T), mean(dat2$resp[dat2$t1==0],na.rm=T),
mean(dat2$viva[dat2$t1==0],na.rm=T), mean(dat2$vic[dat2$t1==0],na.rm=T))
tab <- data.frame('means0' = means0, 'means1' = means1, 'diff' = ates/means0)
rownames(tab) <- c('Messages', 'Respondents', 'Slogans', 'Symbols')
xtable(tab)
#####
## Figure 4 - ATEs over time
#####
## make an index for each 3-hour window
for (i in seq(3,21,3)){
i2 <- (dat2[,paste0('msgs',i)] - mean(dat2[,paste0('msgs',i)], na.rm=T))/sd(dat2[,paste0('msgs',i)][dat2$t1==0], na.rm=T)
i3 <- (dat2[,paste0('resp',i)] - mean(dat2[,paste0('resp',i)], na.rm=T))/sd(dat2[,paste0('resp',i)][dat2$t1==0], na.rm=T)
i4 <- (dat2[,paste0('viva',i)] - mean(dat2[,paste0('viva',i)], na.rm=T))/sd(dat2[,paste0('viva',i)][dat2$t1==0], na.rm=T)
i5 <- (dat2[,paste0('vic',i)] - mean(dat2[,paste0('vic',i)], na.rm=T))/sd(dat2[,paste0('vic',i)][dat2$t1==0], na.rm=T)
dat2$temp <- (i2 + i3 + i4 + i5)/5
setnames(dat2, 'temp', paste0('index',i))
i2 <- (dat2[,paste0('msgs',i,'.pre')] - mean(dat2[,paste0('msgs',i,'.pre')], na.rm=T))/sd(dat2[,paste0('msgs',i,'.pre')][dat2$t1==0], na.rm=T)
i3 <- (dat2[,paste0('resp',i,'.pre')] - mean(dat2[,paste0('resp',i,'.pre')], na.rm=T))/sd(dat2[,paste0('resp',i,'.pre')][dat2$t1==0], na.rm=T)
i4 <- (dat2[,paste0('viva',i,'.pre')] - mean(dat2[,paste0('viva',i,'.pre')], na.rm=T))/sd(dat2[,paste0('viva',i,'.pre')][dat2$t1==0], na.rm=T)
i5 <- (dat2[,paste0('vic',i,'.pre')] - mean(dat2[,paste0('vic',i,'.pre')], na.rm=T))/sd(dat2[,paste0('vic',i,'.pre')][dat2$t1==0], na.rm=T)
i2 <- ifelse(i2==Inf | i2==-Inf, NA, i2)
i3 <- ifelse(i3==Inf | i3==-Inf, NA, i3)
i4 <- ifelse(i4==Inf | i4==-Inf, NA, i4)
i5 <- ifelse(i5==Inf|-Inf, NA, i5)
dat2$temp <- rowMeans(cbind(i2,i3,i4,i5),na.rm=T)
setnames(dat2, 'temp', paste0('index',i,'.pre'))
}
## calculate ITTs
dat2$index24 <- dat2$index.m
outs = c(paste0('index', seq(3,24,3)))
treats = c('t1')
tab <- expand.grid('outs' = outs, 'treats' = treats)
tab$outs <- as.character(tab$outs)
tab$ate <- NA
tab$p <- NA
tab$se_lo <- NA
tab$se_hi <- NA
for (i in 1:dim(tab)[1]){
dat2$outcome <- dat2[,tab$outs[i]]
mod <- lm(as.formula(paste0(tab$outs[i], ' ~ t1 + block + members_log + any_bin + wgt_z')),
data = dat2, weights = ipw)
tab$ate[i] <- summary(mod)$coefficients['t1',1]
tab$p[i] <- summary(mod)$coefficients['t1',4]
tab$se_lo[i] <- summary(mod)$coefficients['t1',1]-1.96*summary(mod)$coefficients['t1',2]
tab$se_hi[i] <- summary(mod)$coefficients['t1',1]+1.96*summary(mod)$coefficients['t1',2]
}
tab$period <- seq(3, 24, 3)
pdf('03_graphs/ate_time_index.pdf')
ggplot(tab, aes(period, ate)) +
geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
geom_errorbar(aes(x = period, y = ate, ymin = se_lo, ymax = se_hi),
lwd = 1, position = position_dodge(width = 1/2),
width = 0.35, colour = "red") +
geom_point(aes(x = period, y = ate),
position = position_dodge(width=1/2), colour = "red") +
theme_minimal() +
theme(text = element_text(size = 20),
axis.text.x = element_text(angle = 0, hjust = 0.5)) +
scale_y_continuous(name = 'Standardized Coefficient') +
scale_x_continuous(name = "Hours Post-Treatment", seq(3, 24, 3))
dev.off()
ggplot(tab, aes(period, ate)) +
geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
geom_errorbar(aes(x = period, y = ate, ymin = se_lo, ymax = se_hi),
lwd = 1, position = position_dodge(width = 1/2),
width = 0.35, colour = "red") +
geom_point(aes(x = period, y = ate),
position = position_dodge(width=1/2), colour = "red") +
theme_minimal() +
theme(text = element_text(size = 20),
axis.text.x = element_text(angle = 0, hjust = 0.5)) +
scale_y_continuous(name = 'Standardized Coefficient') +
scale_x_continuous(name = "Hours Post-Treatment", seq(3, 24, 3))
#####
## Table A3 (appendix) - Heterogeneous effects
#####
dat2$members_log_st <- st(dat2$members_log); dat2$membersXt1 <- dat2$members_log_st * dat2$t1
## models
dat2$pre <- st(dat2$index.pre.m)
mod0a <- lm(index.m ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block,
data = dat2, weights = ipw)
mod0a$se <- coeftest(mod0a, vcov = vcovHC, type = "HC1")[,2]
mod0b <- lm(index.m ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block
+ members_log_st + pre, data = dat2, weights = ipw)
mod0b$se <- coeftest(mod0b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$msgs.pre
mod1a <- lm(msgs ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block,
data = dat2, weights = ipw)
mod1a$se <- coeftest(mod1a, vcov = vcovHC, type = "HC1")[,2]
mod1b <- lm(msgs ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block
+ members_log_st + pre, data = dat2, weights = ipw)
mod1b$se <- coeftest(mod1b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$resp.pre
mod2a <- lm(resp ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block,
data = dat2, weights = ipw)
mod2a$se <- coeftest(mod2a, vcov = vcovHC, type = "HC1")[,2]
mod2b <- lm(resp ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block
+ members_log_st + pre, data = dat2, weights = ipw)
mod2b$se <- coeftest(mod2b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$viva.pre
mod3a <- lm(viva ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block,
data = dat2, weights = ipw)
mod3a$se <- coeftest(mod3a, vcov = vcovHC, type = "HC1")[,2]
mod3b <- lm(viva ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block
+ members_log_st + pre, data = dat2, weights = ipw)
mod3b$se <- coeftest(mod3b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$vic.pre
mod4a <- lm(vic ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block,
data = dat2, weights = ipw)
mod4a$se <- coeftest(mod4a, vcov = vcovHC, type = "HC1")[,2]
mod4b <- lm(vic ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block
+ members_log_st + pre, data = dat2, weights = ipw)
mod4b$se <- coeftest(mod4b, vcov = vcovHC, type = "HC1")[,2]
stargazer(mod0a, mod0b, mod1a, mod1b, mod2a, mod2b, mod3a, mod3b, mod4a, mod4b,
se = list(mod0a$se, mod0b$se, mod1a$se, mod1b$se, mod2a$se, mod2b$se, mod3a$se, mod3b$se, mod4a$se, mod4b$se),
omit = c("block"),
covariate.labels = c('Anger Appeal', 'Power Treatment', 'Any Violence',
'Poverty', 'Members (Log)', 'Pre-Treat Outcome',
'Anger Appeal $\\times$ Power Treatment',
'Anger Appeal $\\times$ Poverty',
'Anger Appeal $\\times$ Any Violence',
'Constant'),
no.space = T, digits = 2,
keep.stat = c('n', 'rsq'))
#####
## Figure 5 - Marginal effects plots
#####
dat2$viol <- dat2$any_bin; dat2$violXt1 <- dat2$viol * dat2$t1
dat2$poorXt1 <- dat2$poor * dat2$t1
dat2$members_log_st <- st(dat2$members_log); dat2$membersXt1 <- dat2$members_log_st * dat2$t1
dat2$t1Xt2 <- dat2$t1 * dat2$t2
dat2$pre <- dat2$index.pre.m
pdf('03_graphs/marg_viol_index.pdf')
print(inter.binning(Y = 'index.m', D = "t1", X = 'viol',
data = dat2,  cutoffs = c(0.5, 1.5),
Z = c('poor', 't2', 't1Xt2', 'poorXt1', 'pre', 'block'),
weights = 'ipw',
na.rm = T, vartype = 'robust', figure = T,
Ylabel = "Mean Effects Index", Dlabel = "Anger Treatment", Xlabel="Violence",
ylab = "Marginal Effect of Anger Treatment",
show.grid = T, theme.bw = T, cex.lab = 1.5, cex.axis = 1.5))
dev.off()
print(inter.binning(Y = 'index.m', D = "t1", X = 'viol',
data = dat2,  cutoffs = c(0.5, 1.5),
Z = c('poor', 't2', 't1Xt2', 'poorXt1', 'pre', 'block'),
weights = 'ipw',
na.rm = T, vartype = 'robust', figure = T,
Ylabel = "Mean Effects Index", Dlabel = "Anger Treatment", Xlabel="Violence",
ylab = "Marginal Effect of Anger Treatment",
show.grid = T, theme.bw = T, cex.lab = 1.5, cex.axis = 1.5))
pdf('03_graphs/marg_poor_index.pdf')
print(inter.binning(Y = 'index.m', D = "t1", X = 'poor',
data = dat2, cutoffs = c(0,1),
Z = c('viol', 't2', 't1Xt2', 'violXt1', 'pre', 'block'),
weights = 'ipw',
na.rm = T, vartype = 'robust', figure = T,
Ylabel = "Mean Effects Index", Dlabel = "Anger Treatment", Xlabel="Poverty",
ylab = "Marginal Effect of Anger Treatment",
show.grid = T, theme.bw = T, cex.lab = 1.5, cex.axis = 1.5))
dev.off()
pdf('03_graphs/marg_t2_index.pdf')
print(inter.binning(Y = 'index.m', D = "t1", X = 't2',
data = dat2, cutoffs = c(0.5, 1.5),
Z = c('viol', 'violXt1', 'poor', 'poorXt1', 'pre', 'block'),
weights = 'ipw',
na.rm = T, vartype = 'robust', figure = T,
Ylabel = "Mean Effects Index", Dlabel = "Anger Treatment", Xlabel="Power Treatment",
ylab = "Marginal Effect of Anger Treatment",
show.grid = T, theme.bw = T, cex.lab = 1.5, cex.axis = 1.5))
dev.off()
#####
## Figure 1: Map constituencies of groups
#####
## read in shape file
shp08 <- readShapePoly('01_data/03_maps/Constituencies.shp') # to open shape file
## export data table
df <- shp08@data
df$gis.id <- 1:nrow(df)
## make list of constits in experiments
sample <- subset(dat2, select = c(constit, province))
sample$harare <- ifelse(dat2$province=="Harare", 1, 0)
sample$constit <- gsub(' [0-9]$', "", sample$constit)
sample$constit <- toupper(sample$constit)
## change names of constits
sample$constit <- gsub(" - ", " ", sample$constit)
df$constit <- gsub(".", "", df$constit, fixed = T)
df$constit <- gsub("'", "", df$constit, fixed = T)
df$constit <- gsub("-", " ", df$constit, fixed = T)
sample$constit <- gsub('ZIBAGWE', 'CHIRUMANZU ZIBAGWE', sample$constit)
sample$constit <- gsub('MUREHWA NORTH', 'MUREWA NORTH', sample$constit)
sample$constit <- gsub('MPOPOMA', 'PELANDABA MPOPOMA', sample$constit)
sample$constit <- gsub('PLUMTREE', 'BEITBRIDGE EAST', sample$constit)
sample$constit <- gsub('WATERFALLS', 'HATFIELD', sample$constit)
## subset to unique
sample <- unique(subset(sample, select = c('constit', 'harare'), is.na(sample$harare)==F))
sample$constit[!(sample$constit %in% df$constit)]
df$constit[!(df$constit %in% sample$constit)]
## make sample indicator
df$sample <- ifelse(df$constit %in% sample$constit, 1, 0)
## make harare indicator
harare <- c("BUDIRIRO", "CHITUNGWIZA SOUTH", "CHITUNGWIZA NORTH", "EPWORTH", "DZIVARASEKWA",
"GLEN NORAH", "GLENVIEW NORTH", "GLENVIEW SOUTH", "HARARE CENTRAL",
"HARARE EAST", "HARARE NORTH", "HARARE SOUTH", "HARARE WEST", "HATFIELD",
"HIGHFIELD EAST", "HIGHFIELD WEST", "KAMBUZUMA", "KUWADZANA",
"KUWADZANA EAST", "MABVUKU TAFARA", "MBARE", "MOUNT PLEASANT",
"MUFAKOSE", "SOUTHERTON", "ST MARYS", "SUNNINGDALE", "WARREN PARK",
"ZENGEZA EAST", "ZENGEZA WEST")
df$harare <- ifelse(df$constit %in% harare, 1, 0)
## make bulawayo indicator
byo <- c("BULAWAYO CENTRAL", "BULAWAYO EAST", "BULAWAYO SOUTH", "EMAKHANDENI-ENTUMBANE",
"LOBENGULA", "LUVEVE", "MAGWEGWE", "MAKOKOBA", "NKETA", "NKULUMANE",
"PELANDABA-MPOPOMA", "PUMULA")
df$byo <- ifelse(df$constit %in% byo, 1, 0)
## put back in shape
df <- df[order(df$gis.id),]
shp08@data <- df
## add province borders to map
prov <- readShapePoly("01_data/03_maps/zim_provinces.shp")
## plot whole zimbabwe
plot_zim <- ggplot() +
layer_spatial(shp08, aes(fill = factor(sample)), show.legend = FALSE) +
scale_fill_manual(values = alpha(c("white", "grey"))) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank()) +
layer_spatial(prov, lwd = 1.2, fill = NA)
## plot harare and bulawayo as insets
har08 <- subset(shp08, shp08$harare==1)
plot_har <- ggplot() +
annotation_spatial(har08) +
layer_spatial(har08, aes(fill = factor(sample)), show.legend = FALSE) +
scale_fill_manual(values = alpha(c("white", "grey")))+
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
byo08 <- subset(shp08, shp08$byo==1)
plot_bul <- ggplot() +
annotation_spatial(byo08) +
layer_spatial(byo08, aes(fill = factor(sample)), show.legend = FALSE) +
scale_fill_manual(values = alpha(c("white", "grey")))+
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
pdf('03_graphs/map_sample_all.pdf')
ggdraw() +
draw_plot(plot_zim) +
theme(plot.margin = unit(c(0,4.5,0,0), "cm")) +
draw_plot(plot_har, x = 1, y = 0.48, width = .3, height = .3) +
draw_plot(plot_bul, x = 1, y = 0.15, width = .3, height = .3) +
geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.5, ymax = 0.8), fill = NA, col = 'black') +
geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.2, ymax = 0.45), fill = NA, col = 'black') +
geom_text(aes(x = 1.15, y = 0.77, label = "Harare"), size = 4) +
geom_text(aes(x = 1.15, y = 0.42, label = "Bulawayo"), size = 4)
dev.off()
#####
## Figure 3: Maps of violence by constituency
#####
sok <- read.csv('01_data/01_raw/sokwanele_violence.csv')
## change sok constit names
setnames(sok, 'consitname', 'constit')
sok$constit <- toupper(sok$constit)
sok$constit <- gsub("MT ", "MOUNT ", sok$constit)
sok$constit <- gsub("GLEN VIEW", "GLENVIEW", sok$constit)
sok$constit <- gsub(" - ", "-", sok$constit)
sok$constit <- gsub("MARAMBA-PFUNGWE", "MARAMBA PFUNGWE", sok$constit)
sok$constit <- gsub("CHIRUMANZU-ZIBAGWE", "CHIRUMANZU ZIBAGWE", sok$constit)
sok$constit <- gsub("ST MARYS", "ST. MARY'S", sok$constit)
sok$constit <- gsub("MABVUKU/TAFARA", "MABVUKU-TAFARA", sok$constit)
sok$constit <- gsub("GOKWE SENGWA", "GOKWE-SENGWA", sok$constit)
sok$constit <- gsub("GOKWE SESAME", "GOKWE-SESAME", sok$constit)
## subset to unique province-constit
constits <- unique(subset(sok, select = c('province', 'constit')))
dup <- duplicated(constits$constit)
constits <- subset(constits, dup==F)
sok$province=NULL
sok <- merge(sok, constits, by = 'constit')
## collapse sok to constit
viol_list <- c('assaultedwithhands', 'assaultedwithweapon',
'abducted', 'falanga', 'submerged', 'strangled', 'burnt',
'hung', 'tortured',
'theftdestrctionproperty')
sok$events <- rowSums(sok[ ,viol_list])
sok$events <- ifelse(sok$events>0,1,0)
sok <- data.table(sok)
sok <- sok[,list('cio'=sum(cid, cio), 'zrp'=sum(zrpuniformed, zrpplainclothes, zrpriot),
'warvet' = sum(zanupfwarveteran), 'supporter' = sum(zanupfsupporter),
'youth' = sum(zanupfyouth), 'zna'=sum(zna),
'events' = sum(events)),
by = c('constit', 'province')]
## read in shape file
zim <- readShapePoly('01_data/03_maps/Constituencies.shp')
## export data table
df <- zim@data
df$gis.id <- 1:nrow(df)
## merge w sok data
df$constit[!(df$constit %in% sok$constit)]
sok$constit[!(sok$constit %in% df$constit)]
df <- merge(df, sok, by = 'constit', all.x = T)
df <- df[order(df$gis.id),]
df$cio <- ifelse(is.na(df$cio)==T, 0, df$cio)
df$zrp <- ifelse(is.na(df$zrp)==T, 0, df$zrp)
df$warvet <- ifelse(is.na(df$warvet)==T, 0, df$warvet)
df$supporter <- ifelse(is.na(df$supporter)==T, 0, df$supporter)
df$youth <- ifelse(is.na(df$youth)==T, 0, df$youth)
df$zna <- ifelse(is.na(df$zna)==T, 0, df$zna)
## calculate proportion
df$all <- df$cio + df$zrp + df$warvet + df$supporter + df$youth + df$zna
df$cio.prop <- df$cio/df$all
df$zrp.prop <- df$zrp/df$all
df$warvet.prop <- df$warvet/df$all
df$supporter.prop <- df$supporter/df$all
df$youth.prop <- df$youth/df$all
df$zna.prop <- df$zna/df$all
## create harare and bulawayo dummies
df$hre <- ifelse(df$constit %in% harare, 1, 0)
df$byo <- ifelse(df$constit %in% byo, 1, 0)
## set violent events to 0 if missing
df$events[is.na(df$events)] <- 0
## put back in map
zim@data <- df
## set colors
brks <- quantile(zim$events, seq(0,1,0.1), na.rm=T)
zim$events_cut <- as.numeric(cut(zim$events, breaks = 10))
cols <- c('white', rev(gray.colors(15)[1:length(brks)]))
## create province layer
prov.layer <- list("sp.lines", prov, col = "black", lwd=1.2)
## plot whole zimbabwe
plot_zim <- ggplot() +
layer_spatial(zim, aes(fill = events)) +
scale_fill_gradientn(colors = cols,
guide = guide_colourbar(title.position = 'top', hjust = 0.5,
barwidth = 15)) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = 'bottom',
legend.text = element_text(size = 14),
legend.title = element_text(size = 16)) +
# scale_shape(guide = guide_legend(title.position = 'top')) +
labs(fill = "Violent Events") +
layer_spatial(prov, lwd = 1.1, fill = NA)
## plot harare and bulawayo as insets
har08 <- subset(zim, zim$hre==1)
plot_har <- ggplot() +
layer_spatial(har08, aes(fill = events_cut), show.legend = FALSE) +
scale_fill_gradientn(colors = cols[1:4]) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
byo08 <- subset(zim, zim$byo==1)
plot_bul <- ggplot() +
layer_spatial(byo08, aes(fill = events_cut), show.legend = FALSE) +
scale_fill_gradientn(colors = cols[1:2]) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
pdf('03_graphs/viol_map.pdf')
ggdraw() +
draw_plot(plot_zim) +
theme(plot.margin = unit(c(0,4.5,0,0), "cm")) +
draw_plot(plot_har, x = 1, y = 0.48, width = .3, height = .3) +
draw_plot(plot_bul, x = 1, y = 0.15, width = .3, height = .3) +
geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.5, ymax = 0.8), fill = NA, col = 'black') +
geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.2, ymax = 0.45), fill = NA, col = 'black') +
geom_text(aes(x = 1.15, y = 0.77, label = "Harare"), size = 5) +
geom_text(aes(x = 1.15, y = 0.42, label = "Bulawayo"), size = 5)
dev.off()
#####
## Figure 3: maps of poverty by constituency
#####
dhs <- read.csv('01_data/02_clean/dhs_clean.csv')
dhs$constit2 <- toupper(dhs$constit2)
dhs$constit2 <- gsub('EMAKHANDENI ENTUMBANE', 'EMAKHANDENI-ENTUMBANE', dhs$constit2)
dhs$constit2 <- gsub('GOKWE ', 'GOKWE-', dhs$constit2)
dhs$constit2 <- gsub('MHONDORO ', 'MHONDORO-', dhs$constit2)
dhs$constit2 <- gsub('ZVISHAVANE ', 'ZVISHAVANE-', dhs$constit2)
dhs$constit2 <- gsub('MUREHWA', 'MUREWA', dhs$constit2)
dhs$constit2 <- gsub('ST MARYS', "ST. MARY'S", dhs$constit2)
dhs$constit2 <- gsub('PELANDABA MPOPOMA', 'PELANDABA-MPOPOMA', dhs$constit2)
dhs$constit2 <- gsub('MABVUKU TAFARA', 'MABVUKU-TAFARA', dhs$constit2)
dhs$constit2[!(dhs$constit2 %in% df$constit)]
df$constit[!(df$constit %in% dhs$constit2)]
df <- merge(df, dhs, by.x = 'constit', by.y = 'constit2', all.x = T)
df$wgt_cuts <- as.numeric(cut(df$wgt_z/100, breaks = 10, include.lowest = T))
cols <- gray.colors(10)
df <- df[order(df$gis.id),]
zim@data <- df
## plot whole zimbabwe
plot_zim <- ggplot() +
layer_spatial(zim, aes(fill = wgt_z/100)) +
scale_fill_gradientn(colors = cols,
guide = guide_colourbar(title.position = 'top', hjust = 0.5,
barwidth = 15)) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = 'bottom',
legend.text = element_text(size = 14),
legend.title = element_text(size = 16)) +
labs(fill = "Weight-for-height z-score") +
layer_spatial(prov, lwd = 1.1, fill = NA)
## plot harare and bulawayo as insets
har08 <- subset(zim, zim$hre==1)
plot_har <- ggplot() +
layer_spatial(har08, aes(fill = wgt_cuts), show.legend = FALSE) +
scale_fill_gradientn(colors = cols[1:9]) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
byo08 <- subset(zim, zim$byo==1)
plot_bul <- ggplot() +
layer_spatial(byo08, aes(fill = wgt_cuts), show.legend = FALSE) +
scale_fill_gradientn(colors = cols[c(4,5,7,10)]) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank())
pdf('03_graphs/pov_map.pdf')
ggdraw() +
draw_plot(plot_zim) +
theme(plot.margin = unit(c(0,4.5,0,0), "cm")) +
draw_plot(plot_har, x = 1, y = 0.48, width = .3, height = .3) +
draw_plot(plot_bul, x = 1, y = 0.15, width = .3, height = .3) +
geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.5, ymax = 0.8), fill = NA, col = 'black') +
geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.2, ymax = 0.45), fill = NA, col = 'black') +
geom_text(aes(x = 1.15, y = 0.77, label = "Harare"), cex = 5) +
geom_text(aes(x = 1.15, y = 0.42, label = "Bulawayo"), cex = 5)
dev.off()
